
# set working directory
setwd("~/Research/Mycorrhizae/Fieldwork/Isotopes")

library(lme4)
library(nlme)
library(lmerTest)
library(olsrr)
library(MuMIn)
library(htmlTable)
library(ggplot2)
library(Rmisc)
library(SciViews)
library(MASS)
library(gridExtra)
library(tidyr)
library(regclass)

# define predictive R^2 functions for later
PRESS <- function(linear.model) {
  #' calculate the predictive residuals
  pr <- residuals(linear.model)/(1-lm.influence(linear.model)$hat)
  #' calculate the PRESS
  PRESS <- sum(pr^2)
  
  return(PRESS)
}

pred_r_squared <- function(linear.model) {
  #' Use anova() to get the sum of squares for the linear model
  lm.anova <- anova(linear.model)
  #' Calculate the total sum of squares
  tss <- sum(lm.anova$'Sum Sq')
  # Calculate the predictive R^2
  pred.r.squared <- 1-PRESS(linear.model)/(tss)
  
  return(pred.r.squared)
}

#########################################
# CREATING GRAPHS OF RAW DATA AND MEANS
################################################

plant.data <- read.csv("WSU plant isotopes for analysis.csv")

# make distance numeric
plant.data$Distance <- as.numeric(plant.data$Distance)

# put distance 5 into 6 bin for analysis
for (i in 1:length(plant.data$Distance)){
  if (plant.data$Distance[i]==5){
    plant.data$Distance[i] <- 6
  }
}

# by stream 
plant.data <- plant.data[ which(plant.data$Stream=="Hansen"),]
plant.data <- plant.data[ which(plant.data$Stream=="Happy"),]
plant.data <- plant.data[ which(plant.data$Stream=="Yako"),]

# by hansen side
plant.data <- plant.data[ which(plant.data$Side=="SL"),]
plant.data <- plant.data[ which(plant.data$Side=="SR"),]

# white spruce or paper birch
plant.data <- plant.data[ which(plant.data$Species=="White Spruce"),]
plant.data <- plant.data[ which(plant.data$Species=="Paper Birch"),]

# choose distance
plant.data <- plant.data[ which(plant.data$Distance < 7),]

# choose transect
plant.data <- plant.data[ which(plant.data$Transect=="S6"),]

# examine by Hansen bank 
tgc <- summarySE(plant.data, measurevar="dC13", groupvars=c("Side"))
tgc
ggplot(tgc, aes(x=Distance, y=total.N, color=Species)) + 
  geom_point(position = position_dodge(.9), stat = "summary", fun = "mean") + 
  geom_line() +
  geom_errorbar(aes(ymin=total.N-se, ymax=total.N+se), width=.2,position=position_dodge(.9)) +
  theme_classic() + 
  geom_point(data=plant.data[which(plant.data$Side=="SL"),], aes(x=Distance, y=total.N), 
             alpha=0.2)  +ylim(0.0091, 0.0274) + scale_x_reverse()

# mean and errors
p <- ggplot(tgc, aes(x=Side, y=dC13)) + 
  geom_errorbar(aes(ymin=dC13-se, ymax=dC13+se), width=.1) + geom_line() + 
  geom_point() + theme_bw() + scale_x_discrete(labels=c('Enhanced', 'Depleted')) + 
  ylab("White spruce dC13") 
p
ggsave(plot=p, width=3, height=3, dpi=300, filename = "dC13_Hansen_spruce_byside.jpg")


# Hansen
# for C:N, use ylim(15,69) 
# for d15N, use ylim(-5,11)
# for d13C, use ylim(-33.5,-28.5)
# for total N, use ylim(0.0091, 0.0274)

# examine at Happy
tgc <- summarySE(plant.data, measurevar="C.N", groupvars=c("Distance","Species"))
ggplot(tgc, aes(x=Distance, y=C.N, color=Species)) + geom_point(position = position_dodge(.9), 
   stat = "summary", fun = "mean") + geom_line() +
  geom_errorbar(aes(ymin=C.N-se, ymax=C.N+se), width=.2,position=position_dodge(.9)) +
  theme_classic() + 
  geom_point(data=plant.data, aes(x=Distance, y=C.N), alpha=0.2)  +
  ylim(16.5,56)

# Happy
# for C:N, use ylim(16.5,56)
# for d15N, use ylim(-3,3.5)
# for dC13, use ylim(-33,-28)

# Yako
# for C:N, use ylim(17,55)
# for d15N, use ylim(-4,6)
#######################################################################
# Plant C:N, dN15 and dC13 by distance and bank at Hansen
#####################################################################
plant.data <- read.csv("WSU plant isotopes for analysis.csv")

# make distance numeric
plant.data$Distance <- as.numeric(plant.data$Distance)

# by stream 
plant.data <- plant.data[ which(plant.data$Stream=="Hansen"),]

# by plant 
plant.data <- plant.data[ which(plant.data$Species=="White Spruce"),]
plant.data <- plant.data[ which(plant.data$Species=="Paper Birch"),]

# select only 1-20 m
plant.data <- plant.data[ which(plant.data$Distance < 21),]

# by plant 
plant.data.ws <- plant.data[ which(plant.data$Species=="White Spruce"),]
plant.data.pb <- plant.data[ which(plant.data$Species=="Paper Birch"),]

# is distance normally distributed

# center Distance variable
#plant.data$Distance <- scale(plant.data$Distance, scale=FALSE)

# create ln of Distance column
plant.data$lnDistance <- ln(plant.data$Distance)

# create quadratic lnDistance column
#plant.data$lnDistance2 <- plant.data$lnDistance^2
#plant.data$Distance2 <- (plant.data$Distance)^2

##############################################
# C:N #
#########################################
# select data
plant.data <- read.csv("WSU plant isotopes for analysis.csv")

# make distance numeric
plant.data$Distance <- as.numeric(plant.data$Distance)

# by stream 
#plant.data <- plant.data[ which(plant.data$Stream=="Hansen"),]

# by plant 
plant.data <- plant.data[ which(plant.data$Species=="White Spruce"),]
plant.data <- plant.data[ which(plant.data$Species=="Paper Birch"),]

# scale variables
test <- scale(plant.data$organic.GWC, scale=FALSE)
plant.data$organic.GWC <- test[,1]

test1 <- scale(plant.data$mineral.GWC, scale=FALSE)
plant.data$mineral.GWC <- test1[,1]

plant.data$decomposing.carcass <- as.factor(plant.data$decomposing.carcass)
plant.data$carcass.deposition <- as.factor(plant.data$carcass.deposition)
plant.data$carcass.load <- as.numeric(plant.data$carcass.load)

# if adding mean mycorrhizal or mean fungal d15N, need to drop rows that do not contain values for these
plant.data <- plant.data %>% drop_na(mean.mycorrhizal.fungal.C.N)
plant.data <- plant.data %>% drop_na(mean.fungal.C.N)

# examine effect of distance on dN15 ratio
# change na. action
options(na.action = "na.fail")
all.parms <- lm(log(C.N) ~ #Side + 
                      ln(Distance) + 
                      #Side:ln(Distance) +
                      #Side:I(ln(Distance)^2) + 
                      #Side:I(ln(Distance)^3) + 
                      #Side:I(ln(Distance)^4) + 
                      mean.fungal.C.N +
                      #mean.mycorrhizal.fungal.C.N +
                      organic.GWC + 
                      mineral.GWC + 
                      decomposing.carcass + 
                      carcass.load + 
                      carcass.deposition, 
                      data=plant.data)
results<-dredge(all.parms)
subset(results, delta <2)
subset(results, delta == 0)

# step model
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

# all data
dredge.model <- lmer(log(C.N) ~ ln(Distance) + organic.GWC + decomposing.carcass + carcass.deposition + carcass.load + 
                       mean.fungal.C.N + (1|Stream) + (1|Transect), data=plant.data)
summary(dredge.model)


dredge.model <- lmer(C.N ~ carcass.load + carcass.deposition + ln(Distance) + mean.fungal.C.N + 
                   mean.mycorrhizal.fungal.C.N + organic.GWC + (1|Stream) + (1|Transect), data=plant.data)
summary(dredge.model)

dredge.model <- lmer(C.N ~ ln(Distance) + (1|Stream:Transect), data=plant.data)
summary(dredge.model)

# get conditional r^2 (which takes into account random effects, second number)
MuMIn::r.squaredGLMM(lme4::lmer(C.N ~ ln(Distance) + organic.GWC + decomposing.carcass + carcass.deposition + carcass.load + 
                                  mean.mycorrhizal.fungal.C.N + (1|Stream), data=plant.data))

# get percentage variance explained by each variable using anova
af <- anova(dredge.model)
afss <- af$"Sum Sq"
blah <- cbind(af,PctExp=afss/sum(afss)*100)
blah$PctExp

# all parameters for a SINGLE SIDE
all.parms<-lm(C.N ~ ln(Distance) +
                I(ln(Distance)^2) + 
                I(ln(Distance)^3) + 
                #I(ln(Distance)^4) + 
                organic.GWC + mineral.GWC,
              data=plant.data)
results<-dredge(all.parms)
subset(results, delta <2)
subset(results, delta == 0)

# get percentage variance explained by each variable using anova
af <- anova(dredge.model)
afss <- af$"Sum Sq"
print(cbind(af,PctExp=afss/sum(afss)*100))

# calculate VIF
library(car)
car::vif(dredge.model)

# 1-20 white spruce C:N with deposition variable
dredge.model <- lm(C.N ~ ln(Distance) + Side + Side:ln(Distance), data=plant.data)
summary(dredge.model)
# decomposing carcass is almost significant

# 1-100 white spruce C:N with deposition variable
dredge.model <- lm(C.N ~ ln(Distance) + Side + carcass.deposition, data=plant.data)
summary(dredge.model)
# decomposing carcass is almost significant

# 1-20 paper birch C:N with deposition variable
dredge.model <- lm(C.N ~ organic.GWC + mineral.GWC + Side, data=plant.data)
summary(dredge.model)
# organic GWC is not significant

# 1-100 paper birch C:N with deposition variable
dredge.model <- lm(C.N ~ ln(Distance) + Side + Side:ln(Distance), data=plant.data)
summary(dredge.model)
# organic GWC is not significant

# for 1-100 white spruce
dredge.model <- lm(C.N ~ Side + Side:I(ln(Distance)^4), data=plant.data)
summary(dredge.model)

# for 1-100 paper birch
dredge.model <- lm(C.N ~ ln(Distance) + Side + Side:ln(Distance) + Side:I(ln(Distance)^4), data=plant.data)
summary(dredge.model)

# for 1-100 paper birch LEFT SIDE
dredge.model <- lm(C.N ~ ln(Distance) + mineral.GWC, data=plant.data)
summary(dredge.model)

# for 1-100 paper birch RIGHT SIDE
dredge.model <- lm(C.N ~ I(ln(Distance)^2) + I(ln(Distance)^4), data=plant.data)
summary(dredge.model)

# for 1-20 white spruce
dredge.model <- lm(C.N ~ decomposing.carcass + mineral.GWC + Side, data=plant.data)
summary(dredge.model)

# for 1-20 paper birch
dredge.model <- lm(C.N ~ mineral.GWC + organic.GWC + Side, data=plant.data)
summary(dredge.model)

# step model
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

pred_r_squared(dredge.model)

# C:N
mod1 <- lm(C.N ~ Side, data=plant.data)
mod2 <- lm(C.N ~ Side + organic.GWC, data=plant.data)
mod3 <- lm(C.N ~ Side + mineral.GWC, data=plant.data)
mod4 <- lm(C.N ~ ln(Distance), data=plant.data)
mod5 <- lm(C.N ~ ln(Distance) + organic.GWC, data=plant.data)
mod6 <- lm(C.N ~ ln(Distance) + mineral.GWC, data=plant.data)
mod7 <- lm(C.N ~ ln(Distance) + Side, data=plant.data)
mod8 <- lm(C.N ~ ln(Distance) + Side + organic.GWC, data=plant.data)
mod9 <- lm(C.N ~ ln(Distance) + Side + mineral.GWC, data=plant.data)
mod10 <- lm(C.N ~ organic.GWC, data=plant.data)
mod11 <- lm(C.N ~ mineral.GWC, data=plant.data)
mod12 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2), data=plant.data)
mod13 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + mineral.GWC, data=plant.data)
mod14 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + organic.GWC, data=plant.data)
mod15 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance), data=plant.data)
mod16 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + mineral.GWC, data=plant.data)
mod17 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + organic.GWC, data=plant.data)
mod18 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3), data=plant.data)
mod19 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.GWC, data=plant.data)
mod20 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + mineral.GWC, data=plant.data)
mod21 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE), data=plant.data)
mod22 <- lm(C.N ~ Side:poly(ln(Distance), 4, raw=TRUE), data=plant.data)
mod23 <- lm(C.N ~ Side:poly(ln(Distance), 5, raw=TRUE), data=plant.data)
mod24 <- lm(C.N ~ Side:poly(ln(Distance), 6, raw=TRUE), data=plant.data)
mod25 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC, data=plant.data)
mod26 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE) + mineral.GWC, data=plant.data)

# for white spruce 1-100, competing top models are step, dredge, and mod21. 
# CHOSE MOD21 - has steep differences between sides and left hump. step and dredge
# are the same model and do not have close bank differences
mod14.lmer <- lmer(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + total.N + (1|Transect), data=plant.data, REML=FALSE)
mod23.lmer <- lmer(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + total.N + organic.GWC + (1|Transect), data=plant.data, REML = FALSE)
mod32.lmer <- lmer(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC + total.N + (1|Transect), data=plant.data, REML=FALSE)

# for white spruce 1-20, all top models have negative predictive r squared, CHOSE STEP
mod5.lmer <- lmer(C.N ~ ln(Distance) + total.N + (1|Transect), data=plant.data, REML=FALSE)
mod6.lmer <- lmer(C.N ~ ln(Distance) + total.N + organic.GWC + (1|Transect), data=plant.data, REML=FALSE)
mod32.lmer <- lmer(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC + total.N + (1|Transect), data=plant.data, REML=FALSE)

# for paper birch 1-100, all models have negative predictive R^2. When modeling
# separately, top model for left side is step.model (positive R^2) and top model
# for right side is dredge but all right models have negative R^2. However, the 
# models for separate sides look basically like the ones for both sides so 
# CHOSE DREDGE 

# for paper birch 1-20, top model is step, CHOSE STEP. all models have negative R^2

results <- AIC(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8, mod9, mod10, mod11, mod12, 
    mod13, mod14, mod15, mod16, mod17, mod18, mod19, mod20, mod21, mod22, mod23, 
    mod24, mod25, mod26, step.model, dredge.model)
#results <- AIC(mod4, mod5, mod6, mod10, mod11, step.model, dredge.model)
ordered <- results[order(-results$AIC),]
ordered
summary(mod17)

select.mod <- step.model

# plot raw data with model estimates
pred_interval <- seq(1, 20, 0.2)
which.side <- "SL"
example <- data.frame(Distance=pred_interval, Side=which.side, 
           total.N = mean(plant.data$total.N[which(plant.data$Side==which.side)]),
           organic.GWC = mean(plant.data$organic.GWC[which(plant.data$Side==which.side)]),
           mineral.GWC = mean(plant.data$mineral.GWC[which(plant.data$Side==which.side)]),
           decomposing.carcass="no")
plant.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
plant.pred1 <- as.data.frame(plant.pred1)
plant.pred1$Distance <- pred_interval

plant.pred1.ws <- plant.pred1
plant.pred1.pb <- plant.pred1

p1 <- ggplot(data=plant.pred1, aes(x=Distance, y=fit))+ geom_line(color="darkgreen")  +
  geom_ribbon(data=plant.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = plant.data[which(plant.data$Side==which.side),], aes(x = Distance, 
  y = C.N, group=Distance, color = NULL), alpha = 0.2, color="limegreen")+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = "black", 
  fill = NA, linewidth = 1.2)) + labs(y = "C.N", x="Distance")+theme_classic()+
  ylim(15.5, 26.5) +
 scale_x_reverse(breaks=c(1,3,6,10,20, 40, 60, 100)) 

# plot raw data with model estimates
pred_interval <- seq(1,20,0.2)
which.side <- "SR"
example <- data.frame(Distance=pred_interval, Side=which.side, 
                      organic.GWC = mean(plant.data$organic.GWC[which(plant.data$Side==which.side)]),
                      mineral.GWC = mean(plant.data$mineral.GWC[which(plant.data$Side==which.side)]),
                      decomposing.carcass="no")
plant.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
plant.pred1 <- as.data.frame(plant.pred1)
plant.pred1$Distance <- pred_interval

plant.pred1.ws <- plant.pred1
plant.pred1.pb <- plant.pred1

p2 <- ggplot(data=plant.pred1, aes(x=Distance, y=fit))+ geom_line(color="darkgreen")  +
  geom_ribbon(data=plant.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = plant.data[which(plant.data$Side==which.side),], aes(x = Distance, 
  y = C.N, group=Distance, color = NULL), alpha = 0.2, color="darkgreen")+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = "black", 
  fill = NA, linewidth = 1.2)) + labs(y = "C.N", x="Distance")+theme_classic()+
  ylim(15.5, 26.5) +
  scale_x_continuous(breaks=c(1,3,6,10,20, 40, 60, 100)) 

grid.arrange(p1, p2, nrow=1)

ggsave(plot=p1, width=5, height=3, dpi=300, filename = "CN_20_paperbirch_L.jpg")
ggsave(plot=p2, width=5, height=3, dpi=300, filename = "CN_20_paperbirch_R.jpg")

# for white spruce 1-100 C:N, use limits ylim(39,69)
# for white spruce 1-20 C:N, use limits ylim(39,57.2)

# for paper birch 1-100 C:N, use limits ylim(15.5, 32)
# for paper birch 1-20 C:N, use limits ylim(15.5, 26.5)

# model diagnostics
# see if residuals are normally distributed
plot(density(mod8$residuals))
qqnorm(mod8$residuals)
qqline(mod8$residuals)
shapiro.test(mod8$residuals)
# plot residuals against predicted values to determine heteroscedacitty
plot(fitted(mod8),mod8$residuals)
abline(0,0)

#####################################################################
# dN15 - here, include total N but do not include GWC
#####################################################################
# select data
plant.data <- read.csv("WSU plant isotopes for analysis.csv")

# make distance numeric
plant.data$Distance <- as.numeric(plant.data$Distance)

# by stream 
plant.data <- plant.data[ which(plant.data$Stream=="Hansen"),]

# by plant 
plant.data <- plant.data[ which(plant.data$Species=="White Spruce"),]
plant.data <- plant.data[ which(plant.data$Species=="Paper Birch"),]

# scale variables
test <- scale(plant.data$X.N, scale=FALSE)
plant.data$X.N <- test[,1]

test1 <- scale(plant.data$organic.GWC, scale=FALSE)
plant.data$organic.GWC <- test1[,1]

test2 <- scale(plant.data$mineral.GWC, scale=FALSE)
plant.data$mineral.GWC <- test2[,1]

test2 <- scale(plant.data$mean.mycorrhizal.fungal.d15N, scale=FALSE)
plant.data$mean.mycorrhizal.fungal.d15N <- test2[,1]

test2 <- scale(plant.data$mean.fungal.d15N, scale=FALSE)
plant.data$mean.fungal.d15N <- test2[,1]

test2 <- scale(plant.data$organic.soil.d15N, scale=FALSE)
plant.data$organic.soil.d15N <- test2[,1]

test2 <- scale(plant.data$mineral.soil.d15N, scale=FALSE)
plant.data$mineral.soil.d15N <- test2[,1]

plant.data$decomposing.carcass <- as.factor(plant.data$decomposing.carcass)
plant.data$carcass.deposition <- as.factor(plant.data$carcass.deposition)
plant.data$carcass.load <- as.numeric(plant.data$carcass.load)

# examine effect of distance on dN15 ratio

# for hump, omit transect and mean fungal d15N

# if adding mean mycorrhizal or mean fungal d15N, need to drop rows that do not contain values for these
plant.data <- plant.data %>% drop_na(mean.mycorrhizal.fungal.d15N)
plant.data <- plant.data %>% drop_na(mean.fungal.d15N)

# for white spruce 1-20, have side and distance variables up to fourth power (or just to second) only
# for paper birch, have squared side and distance, plus X.N

# change na. action
options(na.action = "na.fail")
all.parms<-lm(d15N ~ 
                #Side + 
                ln(Distance) + 
                #Side:ln(Distance) +
                #Side:I(ln(Distance)^2) +  
                #Side:I(ln(Distance)^3) + 
                #Side:I(ln(Distance)^4),
                carcass.deposition + 
                decomposing.carcass + 
                carcass.load + 
                X.N + 
                organic.GWC +
                mineral.GWC + 
                mean.mycorrhizal.fungal.d15N +
                #mean.fungal.d15N +
                #mean.fungal.d15N:Side + 
                organic.soil.d15N + 
                mineral.soil.d15N,
                #mean.mycorrhizal.fungal.d15N:Side + 
                #Transect +
              data=plant.data)
results<-dredge(all.parms)
subset(results, delta <2)
subset(results, delta == 0)

# step model
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

 # step model with random effedcts
step.model <- lmer(d15N ~ carcass.deposition + ln(Distance) + mineral.soil.d15N +
                     mean.fungal.d15N + organic.soil.d15N + 
                     (1|Stream:Transect), data=plant.data)
summary(step.model)

#dredge model
dredge.model <- lmer(d15N ~ carcass.deposition + ln(Distance) + mean.fungal.d15N + mineral.soil.d15N + organic.GWC + organic.soil.d15N +
                       (1|Stream:Transect), data=plant.data)
summary(dredge.model)

# all data model
dredge.model <- lmer(d15N ~  ln(Distance) +carcass.deposition + decomposing.carcass +  carcass.load + X.N + organic.GWC +mineral.GWC + 
                mean.fungal.d15N +organic.soil.d15N + mineral.soil.d15N + (1:Stream|Transect), data=plant.data)
summary(dredge.model)

# model with specific hypotheses
mod1 <- lmer(d15N ~  ln(Distance) +carcass.deposition + decomposing.carcass +  carcass.load + X.N + 
                organic.GWC + organic.soil.d15N + 
                mean.fungal.d15N + (1:Stream|Transect), data=plant.data)
summary(mod1)

mod2 <- lmer(d15N ~  ln(Distance) +carcass.deposition + decomposing.carcass +  carcass.load + X.N + 
                mineral.GWC + mineral.soil.d15N + 
                mean.fungal.d15N + (1:Stream|Transect), data=plant.data)
summary(mod2)

mod3 <- lmer(d15N ~  ln(Distance) +carcass.deposition + decomposing.carcass +  carcass.load + X.N + 
                organic.GWC + organic.soil.d15N + 
                mean.mycorrhizal.fungal.d15N + (1:Stream|Transect), data=plant.data)
summary(mod3)

mod4 <- lmer(d15N ~  ln(Distance) +carcass.deposition + decomposing.carcass +  carcass.load + X.N + 
                mineral.GWC + mineral.soil.d15N + 
                mean.mycorrhizal.fungal.d15N + (1:Stream|Transect), data=plant.data)
summary(mod4)

AIC(mod1, mod2, mod3, mod4)

ggplot(plant.data,aes(x=Transect,y=d15N)) + geom_boxplot(alpha=0.2)

# get conditional r^2 (which takes into account random effects, second number)
MuMIn::r.squaredGLMM(lme4::lmer(d15N ~  ln(Distance) +carcass.deposition + decomposing.carcass +  carcass.load + X.N + organic.GWC +mineral.GWC + 
                                  mean.fungal.d15N +organic.soil.d15N + mineral.soil.d15N + (1:Stream|Transect), data=plant.data))

# paper birch
dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass + ln(Distance) + mean.fungal.d15N + mineral.soil.d15N + 
                       organic.GWC + organic.soil.d15N, data=plant.data)
dredge.model <- lmer(d15N ~ carcass.deposition + decomposing.carcass + ln(Distance) + mean.fungal.d15N + mineral.soil.d15N + 
                       organic.GWC + organic.soil.d15N + (1|Stream:Transect), data=plant.data)
summary(dredge.model)

# get conditional r^2 (which takes into account random effects, second number)
MuMIn::r.squaredGLMM(lme4::lmer(d15N ~ carcass.deposition + decomposing.carcass + ln(Distance) + mean.fungal.d15N + mineral.soil.d15N + 
                                  organic.GWC + organic.soil.d15N + (1|Stream:Transect), data=plant.data))

pred_r_squared(dredge.model)

# get percentage variance explained by each variable using anova
af <- anova(dredge.model)
afss <- af$"Sum Sq"
blah <- cbind(af,PctExp=afss/sum(afss)*100)
blah$PctExp

# calculate VIF
library(car)
car::vif(dredge.model)

# 
dredge.model <- lm(d15N ~ carcass.deposition + organic.GWC + organic.soil.d15N + decomposing.carcass, data=plant.data)
summary(dredge.model)

# white spruce 1-100
dredge.model <- lm(d15N ~ carcass.deposition + mean.fungal.d15N + organic.GWC + organic.soil.d15N, data=plant.data)
summary(dredge.model)

# white spruce 1-20
dredge.model <- lmer(d15N ~ carcass.deposition + mean.fungal.d15N + organic.soil.d15N + (1|Transect), data=plant.data)
summary(dredge.model)

# paper birch 1-100
dredge.model <- lmer(d15N ~ carcass.deposition + mean.fungal.d15N + organic.GWC + organic.soil.d15N + (1|Transect), data=plant.data)
summary(dredge.model)

# paper birch 1-20
dredge.model <- lmer(d15N ~ carcass.deposition + organic.GWC + organic.soil.d15N + (1|Transect), data=plant.data)
summary(dredge.model)

# adding distance gives negative effect of distance but doesn't change anything else for white spruce, for paper birch distance has no effect
dredge.model <- lm(d15N ~ carcass.deposition + ln(Distance) + organic.GWC + organic.soil.d15N, data=plant.data)
summary(dredge.model)

# for white spruce 1-20, top model from dredge , weak association with MF (p=0.08) (side and distance not included, when included, it doesn't show MF or organic soil as result)
dredge.model <- lm(d15N ~ carcass.deposition + mean.mycorrhizal.fungal.d15N + organic.soil.d15N, data=plant.data)
summary(dredge.model)

# for white spruce 1-100, top model from dredge
dredge.model <- lm(d15N ~ carcass.deposition + X.N + organic.GWC + organic.soil.d15N, data=plant.data)
summary(dredge.model)

# for paper birch 1-20, top model from dredge (side and distance included)
dredge.model <- lm(d15N ~ carcass.deposition + organic.GWC + organic.soil.d15N, data=plant.data)
summary(dredge.model)

# for paper birch 1-100, top model from dredge (side and distance included)
dredge.model <- lm(d15N ~ carcass.deposition + organic.GWC + organic.soil.d15N, data=plant.data)
summary(dredge.model)

# for white spruce 1-100, top model from dredge 
dredge.model <- lm(d15N ~ carcass.deposition + organic.GWC + organic.soil.d15N + X.N, data=plant.data)
# for white spruce 1-100, top model from dredge within 2 AIC and this is also all of the variables
dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass + mineral.soil.d15N + organic.GWC + organic.soil.d15N + X.N, data=plant.data)
# both yielded same significant variables - carcass deposition (+), organic GWC (+), and organic soil d15N (+)

# for white spruce 1-100, top model from dredge with MF
dredge.model <- lm(d15N ~ carcass.deposition + organic.GWC + organic.soil.d15N + X.N, data=plant.data)
# for white spruce 1-100, top model from dredge within 2 AIC 
dredge.model <- lm(d15N ~ carcass.deposition + organic.GWC + organic.soil.d15N + X.N + mineral.soil.d15N + mean.mycorrhizal.fungal.d15N, data=plant.data)
# all variables with MF
dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass + mean.mycorrhizal.fungal.d15N + organic.GWC + organic.soil.d15N + X.N + 
                     mean.mycorrhizal.fungal.d15N:Side + mineral.soil.d15N, data=plant.data)
# all yielded same significant variables - carcass deposition (+), organic GWC (+), and organic soil d15N (+)

# for paper birch 1-100, top model from dredge 
dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass + organic.GWC + organic.soil.d15N, data=plant.data)
# for paper birch 1-100, top model from dredge within 2 AIC
dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass + organic.GWC + organic.soil.d15N + mineral.soil.d15N, data=plant.data)
# all variables
dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass + mineral.soil.d15N + organic.GWC + organic.soil.d15N + X.N, data=plant.data)
# both yielded same significant variables - carcass deposition (+), organic GWC (+), and organic soil d15N (+)
summary(dredge.model)
# all yielded same significant variables - carcass deposition (+), organic GWC (+), and organic soil d15N (+)

# for paper birch 1-100, top model from dredge with MF
dredge.model <- lm(d15N ~ carcass.deposition + organic.GWC + organic.soil.d15N + mean.mycorrhizal.fungal.d15N, data=plant.data)
# for paper birch 1-100, top model from dredge within 2 AIC 
dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass + mean.mycorrhizal.fungal.d15N + organic.GWC + organic.soil.d15N +  
                     mean.mycorrhizal.fungal.d15N:Side, data=plant.data)
# all variables with MF
dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass + mean.mycorrhizal.fungal.d15N + organic.GWC + organic.soil.d15N + X.N + 
                     mean.mycorrhizal.fungal.d15N:Side + mineral.soil.d15N, data=plant.data)
# all yielded same significant variables - carcass deposition (+), organic GWC (+), and organic soil d15N (+)


# for white spruce 1-100 
dredge.model <- lm(d15N ~ decomposing.carcass + Side + ln(Distance) + mean.mycorrhizal.fungal.d15N + X.N + mean.mycorrhizal.fungal.d15N:Side, data=plant.data)
summary(dredge.model)

# for paper birch 1-100
dredge.model <- lm(d15N ~  ln(Distance) + mean.mycorrhizal.fungal.d15N + organic.GWC + Side + Side:ln(Distance) + mean.mycorrhizal.fungal.d15N:Side, data=plant.data)

# for white spruce 1-20
dredge.model <- lm(d15N ~ ln(Distance) + mean.mycorrhizal.fungal.d15N + mean.mycorrhizal.fungal.d15N:Side + Side, data=plant.data)
summary(dredge.model)

# for paper birch 1-20
dredge.model <- lm(d15N ~ Side + Side:I(ln(Distance)^2), data=plant.data)
summary(dredge.model)

step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

pred_r_squared(dredge.model)

# ln distance
mod1 <- lm(d15N ~ Side, data=plant.data)
mod2 <- lm(d15N ~ Side + X.N, data=plant.data)
mod3 <- lm(d15N ~ Side + X.N + organic.GWC, data=plant.data)
mod4 <- lm(d15N ~ Side + X.N + mineral.GWC, data=plant.data)
mod5 <- lm(d15N ~ ln(Distance) + X.N, data=plant.data)
mod6 <- lm(d15N ~ ln(Distance) + X.N + organic.GWC, data=plant.data)
mod7 <- lm(d15N ~ ln(Distance) + X.N + mineral.GWC, data=plant.data)
mod8 <- lm(d15N ~ ln(Distance), data=plant.data)
mod9 <- lm(d15N ~ ln(Distance) + Side, data=plant.data)
mod10 <- lm(d15N ~ ln(Distance) + Side + X.N, data=plant.data)
mod11 <- lm(d15N ~ ln(Distance) + Side + organic.GWC, data=plant.data)
mod12 <- lm(d15N ~ ln(Distance) + Side + mineral.GWC, data=plant.data)
mod13 <- lm(d15N ~ X.N, data=plant.data)
mod14 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N, data=plant.data)
mod15 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + mineral.GWC, data=plant.data)
mod16 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + organic.GWC, data=plant.data)
mod17 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance), data=plant.data)
mod18 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + X.N, data=plant.data)
mod19 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2), data=plant.data)
mod20 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + mineral.GWC, data=plant.data)
mod21 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + organic.GWC, data=plant.data)
mod22 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N + mineral.GWC, data=plant.data)
mod23 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N + organic.GWC, data=plant.data)
mod24 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3), data=plant.data)
mod24a <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.GWC, data=plant.data)
mod24b <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.GWC + X.N, data=plant.data)
mod24c <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + mineral.GWC + X.N, data=plant.data)
mod25 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^3), data=plant.data)
# although model 15s have as low AIC as model 12, they do not include individual
# variables of side and distance which are important 
mod26a <- lm(d15N ~ Side:I(ln(Distance)^2) + X.N, data=plant.data)
mod26 <- lm(d15N ~ Side:I(ln(Distance)^3) + X.N, data=plant.data)
mod27 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE), data=plant.data)
mod28 <- lm(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE), data=plant.data)
mod29 <- lm(d15N ~ Side:poly(ln(Distance), 5, raw=TRUE), data=plant.data)
mod30 <- lm(d15N ~ Side:poly(ln(Distance), 6, raw=TRUE), data=plant.data)
mod31 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC + X.N, data=plant.data)
mod32 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + mineral.GWC + X.N, data=plant.data)
mod33 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + X.N, data=plant.data)
mod34 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC, data=plant.data)
mod35 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + mineral.GWC, data=plant.data)

mod24.lmer <- lmer(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + (1|Transect), data=plant.data, REML=FALSE)

# model comparison
results <- AIC(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8, mod9, mod10, mod11, mod12, mod13,
    mod14, mod15, mod16, mod17, mod18, mod19, mod20, mod21, mod22, mod23, mod24,
    mod24a, mod24b, mod24c, mod25, mod26, mod26a, mod27, mod28, mod29, mod30,
    mod31, mod32, mod33, mod34, mod35, dredge.model, step.model)
ordered <- results[order(-results$AIC),]
ordered
summary(mod18)

# PLANTS COMBINED
# for 1-100, best model is mod24
# for 1-20, best model is mod24 (all other models within 2 AIC are worse)

# PAPER BIRCH
# for 1-100, top model is step
# for 1-20, choose step model

# WHITE SPRUCE
# for 1-100, top is step
# for 1-20, chose dredge (step is overfit)

# fit mean mycorrhizal model
left <- plant.data[which(plant.data$Side=="SL"),]
m.fun <- lm(mean.mycorrhizal.fungal.d15N ~ I(ln(Distance))^2, data=left)
pred_interval <- seq(1,20,0.2)
example <- data.frame(Distance=pred_interval)
m.fun.pred <- predict(m.fun, newdata=example, re.form=NA, type="response", interval='confidence')
plot(m.fun.pred)
m.fun.pred.fit <- m.fun.pred[,1]

c(5.581784, 5.255117, 2.705117, -1.520716)

# select the model that will be used to plot predictions
select.mod <- step.model

# plot raw data with model estimates for LEFT bank
pred_interval <- seq(1,20,0.2)
#pred_interval <- c(1,6,10,20)
which.side <- "SL"
example <- data.frame(Distance=pred_interval, Side=which.side, 
           X.N = mean(plant.data$X.N[which(plant.data$Side==which.side)]),
           organic.GWC = mean(plant.data$organic.GWC[which(plant.data$Side==which.side)]),
           mineral.GWC = mean(plant.data$mineral.GWC[which(plant.data$Side==which.side)]),
           Transect="S6",
           decomposing.carcass="yes",
           mean.mycorrhizal.fungal.d15N = mean(plant.data$mean.mycorrhizal.fungal.d15N[which(plant.data$Side==which.side)]))
           #mean.mycorrhizal.fungal.d15N = m.fun.pred[,1])
plant.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
plant.pred1 <- as.data.frame(plant.pred1)
plant.pred1$Distance <- pred_interval

plant.pred1.ws <- plant.pred1
plant.pred1.pb <- plant.pred1

p1 <- ggplot(data=plant.pred1, aes(x=Distance, y=fit))+ geom_line(color="darkgreen") +
  geom_ribbon(data=plant.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = plant.data[which(plant.data$Side==which.side),], aes(x = Distance, 
  y = d15N, group=Distance), alpha = 0.1, color = "darkgreen")+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = "black", 
  fill = NA, linewidth = 1.2)) + labs(y = "d15N", x="Distance") + 
  theme_classic() + theme(legend.position="none") +
  ylim(-0.5,13) +
  scale_x_reverse(breaks=c(1,3,6,10,20, 40, 60, 100))

# plot raw data with model estimates for RIGHT bank
pred_interval <- seq(1,20,0.2)
which.side <- "SR"
example <- data.frame(Distance=pred_interval, Side=which.side, 
                      X.N = mean(plant.data$X.N[which(plant.data$Side==which.side)]),
                      organic.GWC = mean(plant.data$organic.GWC[which(plant.data$Side==which.side)]),
                      mineral.GWC = mean(plant.data$mineral.GWC[which(plant.data$Side==which.side)]),
                      Transect="S6",
                      decomposing.carcass="yes",
                      mean.mycorrhizal.fungal.d15N = mean(plant.data$mean.mycorrhizal.fungal.d15N[which(plant.data$Side==which.side)]))
plant.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
plant.pred1 <- as.data.frame(plant.pred1)
plant.pred1$Distance <- pred_interval

plant.pred1.ws <- plant.pred1
plant.pred1.pb <- plant.pred1

p2 <- ggplot(data=plant.pred1, aes(x=Distance, y=fit))+ geom_line(color="darkgreen") +
  geom_ribbon(data=plant.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = plant.data[which(plant.data$Side==which.side),], aes(x = Distance, 
  y = d15N, group=Distance), alpha = 0.1, color = "darkgreen")+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = "black", 
  fill = NA, linewidth = 1.2)) + labs(y = "d15N", x="Distance") + 
  theme_classic() + theme(legend.position="none") +
  ylim(-0.5,13) +
scale_x_continuous(breaks=c(1,3,6,10,20, 40, 60, 100))

grid.arrange(p1, p2, nrow=1)

ggsave(plot=p1, width=5, height=3, dpi=300, filename = "dN15_20_paperbirch_L.jpg")
ggsave(plot=p2, width=5, height=3, dpi=300, filename = "dN15_20_paperbirch_R.jpg")

# for plant 1-100, ylim(-4.93, 11.4)
# for plant 1-20, ylim(-1.46, 18.2)

# for white spruce 1-100, ylim(-4.93, 18)
# for white spruce 1-20, ylim(-0.5,15)

# for paper birch 1-100, ylim(-4, 12.5)
# for paper birch 1-20, ylim(-1.5, 18)

# calculate variance inflation factor for correlated variables total N and GWC
library(car)
library(rms)
library(DAAG)
rms::vif(mod12)
car::vif(mod12)
DAAG::vif(mod12)
# looks like it's okay to include both total N and GWC, only .4 correlation anyways

# untransformed distance to compare (can check this later)

##############################################################
# dC13 #
##############################################################
# select data
plant.data <- read.csv("WSU plant isotopes for analysis.csv")

# make distance numeric
plant.data$Distance <- as.numeric(plant.data$Distance)

# by stream 
plant.data <- plant.data[ which(plant.data$Stream=="Hansen"),]

# by plant 
plant.data <- plant.data[ which(plant.data$Species=="White Spruce"),]
plant.data <- plant.data[ which(plant.data$Species=="Paper Birch"),]

# scale variables
test <- scale(plant.data$X.N, scale=FALSE)
plant.data$X.N <- test[,1]

test1 <- scale(plant.data$organic.GWC, scale=FALSE)
plant.data$organic.GWC <- test1[,1]

test2 <- scale(plant.data$mineral.GWC, scale=FALSE)
plant.data$mineral.GWC <- test2[,1]

plant.data$decomposing.carcass <- as.factor(plant.data$decomposing.carcass)
plant.data$carcass.deposition <- as.factor(plant.data$carcass.deposition)
plant.data$carcass.load <- as.numeric(plant.data$carcass.load)

# change na. action
options(na.action = "na.fail")
all.parms<-lm(dC13 ~ 
                #Side + 
                ln(Distance) + 
                #Side:ln(Distance) +
                #Side:I(ln(Distance)^2) + 
                #Side:I(ln(Distance)^3) + 
                #Side:I(ln(Distance)^4) +
                X.N +
                organic.GWC + 
                #mineral.GWC + 
                decomposing.carcass + 
                carcass.deposition,
              data=plant.data)
results<-dredge(all.parms)
subset(results, delta <2)
subset(results, delta == 0)

# all data
dredge.model <- lmer(dC13 ~ ln(Distance) + X.N + organic.GWC + decomposing.carcass + carcass.deposition + 
                       carcass.load + (1|Stream), data=plant.data)
summary(dredge.model)

dredge.model <- lm(dC13 ~ decomposing.carcass + ln(Distance) + X.N, data=plant.data)
summary(dredge.model)

dredge.model <- lmer(dC13 ~ decomposing.carcass + ln(Distance) + X.N + (1|Stream), data=plant.data)
summary(dredge.model)

# get conditional r^2 (which takes into account random effects, second number)
MuMIn::r.squaredGLMM(lme4::lmer(dC13 ~ ln(Distance) + X.N + organic.GWC + decomposing.carcass + carcass.deposition + 
                                  carcass.load + (1|Stream:Transect), data=plant.data))

# get percentage variance explained by each variable using anova
af <- anova(dredge.model)
afss <- af$"Sum Sq"
blah <- cbind(af,PctExp=afss/sum(afss)*100)
blah$PctExp

# calculate VIF
library(car)
car::vif(dredge.model)

# all parameters for a SINGLE SIDE
all.parms<-lm(dC13 ~ ln(Distance) + ln(Distance) +
                I(ln(Distance)^2) + I(ln(Distance)^3) + 
                I(ln(Distance)^4) + X.N + organic.GWC + mineral.GWC,
              data=plant.data)
results<-dredge(all.parms)
subset(results, delta <2)
subset(results, delta == 0)

# get percentage variance explained by each variable using anova
af <- anova(dredge.model)
afss <- af$"Sum Sq"
print(cbind(af,PctExp=afss/sum(afss)*100))

# calculate VIF
library(car)
car::vif(dredge.model)

# for 1-20 white spruce dC13 (with deposition variable)
dredge.model <- lm(dC13 ~ ln(Distance), data=plant.data)
summary(dredge.model)

# for 1-100 white spruce dC13 (with deposition variable)
dredge.model <- lm(dC13 ~ ln(Distance) + decomposing.carcass + Side + X.N, data=plant.data)
summary(dredge.model)

# for 1-20 paper birch dC13 (with deposition variable)
dredge.model <- lm(dC13 ~ ln(Distance) + decomposing.carcass, data=plant.data)
summary(dredge.model)

# for 1-100 paper birch dC13 (with deposition variable)
dredge.model <- lm(dC13 ~ ln(Distance) + decomposing.carcass, data=plant.data)
summary(dredge.model)

# for 1-100 paper birch LEFT AND RIGHT SIDE SEPARATELY
dredge.model <- lm(dC13 ~ 1, data=plant.data)
summary(dredge.model)

# for 1-20 paper birch LEFT AND RIGHT SIDE SEPARATELY
dredge.model <- lm(dC13 ~ 1, data=plant.data)
summary(dredge.model)

# for 1-100 white spruce
dredge.model <- lm(dC13 ~ Side + ln(Distance) + Side:I(ln(Distance)^4), data=plant.data)
summary(dredge.model)

# for 1-20 white spruce
dredge.model <- lm(dC13 ~ ln(Distance), data=plant.data)

# for 1-100 paper birch
dredge.model <- lm(dC13 ~ 1, data=plant.data)
summary(dredge.model)

# for 1-20 paper birch
dredge.model <- lm(dC13 ~ 1, data=plant.data)

step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

pred_r_squared(dredge.model)

# ln distance
mod1 <- lm(dC13 ~ Side, data=plant.data)
mod2 <- lm(dC13 ~ Side + X.N, data=plant.data)
mod3 <- lm(dC13 ~ Side + X.N + organic.GWC, data=plant.data)
mod4 <- lm(dC13 ~ Side + X.N + mineral.GWC, data=plant.data)
mod5 <- lm(dC13 ~ ln(Distance) + X.N, data=plant.data)
mod6 <- lm(dC13 ~ ln(Distance) + X.N + organic.GWC, data=plant.data)
mod7 <- lm(dC13 ~ ln(Distance) + X.N + mineral.GWC, data=plant.data)
mod8 <- lm(dC13 ~ ln(Distance), data=plant.data)
mod9 <- lm(dC13 ~ ln(Distance) + Side, data=plant.data)
mod10 <- lm(dC13 ~ ln(Distance) + Side + X.N, data=plant.data)
mod11 <- lm(dC13 ~ ln(Distance) + Side + organic.GWC, data=plant.data)
mod12 <- lm(dC13 ~ ln(Distance) + Side + mineral.GWC, data=plant.data)
mod13 <- lm(dC13 ~ X.N, data=plant.data)
mod14 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N, data=plant.data)
mod15 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + mineral.GWC, data=plant.data)
mod16 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + organic.GWC, data=plant.data)
mod17 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance), data=plant.data)
mod18 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + X.N, data=plant.data)
mod19 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2), data=plant.data)
mod20 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + mineral.GWC, data=plant.data)
mod21 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + organic.GWC, data=plant.data)
mod22 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N + mineral.GWC, data=plant.data)
mod23 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N + organic.GWC, data=plant.data)
mod24 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3), data=plant.data)
mod24a <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.GWC, data=plant.data)
mod24b <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.GWC + X.N, data=plant.data)
mod24c <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + mineral.GWC + X.N, data=plant.data)
mod25 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^3), data=plant.data)
# although model 15s have as low AIC as model 12, they do not include individual
# variables of side and distance which are important 
mod26a <- lm(dC13 ~ Side:I(ln(Distance)^2) + X.N, data=plant.data)
mod26 <- lm(dC13 ~ Side:I(ln(Distance)^3) + X.N, data=plant.data)
mod27 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE), data=plant.data)
mod28 <- lm(dC13 ~ Side:poly(ln(Distance), 4, raw=TRUE), data=plant.data)
mod29 <- lm(dC13 ~ Side:poly(ln(Distance), 5, raw=TRUE), data=plant.data)
mod30 <- lm(dC13 ~ Side:poly(ln(Distance), 6, raw=TRUE), data=plant.data)
mod31 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC + X.N, data=plant.data)
mod32 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + mineral.GWC + X.N, data=plant.data)
mod33 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + X.N, data=plant.data)
mod34 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC, data=plant.data)
mod35 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + mineral.GWC, data=plant.data)

# model comparison
results <- AIC(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8, mod9, mod10, mod11, mod12, mod13,
               mod14, mod15, mod16, mod17, mod18, mod19, mod20, mod21, mod22, mod23, mod24,
               mod24a, mod24b, mod24c, mod25, mod26, mod26a, mod27, mod28, mod29, mod30,
               mod31, mod32, mod33, mod34, mod35, step.model, dredge.model)
results <- AIC(mod5, mod6, mod7, mod8, mod13, step.model, dredge.model)
ordered <- results[order(-results$AIC),]
ordered
summary(mod18)

# for white spruce 1-100, tons of competing models, can't decide between step mod 
# and mod33, CHOSE MOD33 which is top
# for white spruce 1-20, dredge is top model, CHOSE DREDGE

# for paper birch dC13 1-100, all models have negative predictive R squared 
# when modeling separate banks, left and right bank still had negative R^2. 
# for left bank, best model by far is step, for right best model is mod8. For 1-20m,
# best model is mod8 so maybe should just go with mod8 for 1-100

# for paper birch dC13 1-20, all models have negative predictive R squared and also
# when modeling banks separately so just chose top model which is dredge.model

pred_r_squared(mod8)
# select the model that will be used to plot predictions
select.mod <- dredge.model

# plot raw data with model estimates
pred_interval <- seq(1,20,0.2)
which.side <- "SL"
example <- data.frame(Distance=pred_interval, Side=which.side, 
                      organic.GWC=mean(plant.data$organic.GWC[which(plant.data$Side==which.side)]),
                      mineral.GWC=mean(plant.data$mineral.GWC[which(plant.data$Side==which.side)]),
                      X.N=mean(plant.data$X.N[which(plant.data$Side==which.side)]))
plant.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
plant.pred1 <- as.data.frame(plant.pred1)
plant.pred1$Distance <- pred_interval

plant.pred1.ws <- plant.pred1
plant.pred1.pb <- plant.pred1

p1 <- ggplot(data=plant.pred1, aes(x=Distance, y=fit))+ geom_line(color="darkgreen") + 
  geom_ribbon(data=plant.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = plant.data[which(plant.data$Side==which.side),], aes(x = Distance, 
  y = dC13, group=Distance, color = NULL), alpha = 0.1, color="darkgreen")+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = 
  "black", fill = NA, linewidth = 1.2)) + labs(y = "d13C", x="Distance") + 
  theme_classic() + theme(legend.position="none") +
  ylim(-33.1, -29) + 
  scale_x_reverse(breaks=c(1,3,6,10,20, 40, 60, 100)) 

# plot raw data with model estimates for RIGHT BANK
pred_interval <- seq(1,20,0.2)
which.side <- "SR"
example <- data.frame(Distance=pred_interval, Side=which.side, 
           organic.GWC=mean(plant.data$organic.GWC[which(plant.data$Side==which.side)]),
           mineral.GWC=mean(plant.data$mineral.GWC[which(plant.data$Side==which.side)]),
           X.N=mean(plant.data$X.N[which(plant.data$Side==which.side)]))
plant.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
plant.pred1 <- as.data.frame(plant.pred1)
plant.pred1$Distance <- pred_interval

plant.pred1.ws <- plant.pred1
plant.pred1.pb <- plant.pred1

p2 <- ggplot(data=plant.pred1, aes(x=Distance, y=fit))+ geom_line(color="darkgreen") + 
  geom_ribbon(data=plant.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = plant.data[which(plant.data$Side==which.side),], aes(x = Distance, 
  y = dC13, group=Distance, color = NULL), alpha = 0.1, color="darkgreen")+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = 
  "black", fill = NA, linewidth = 1.2)) + labs(y = "d13C", x="Distance") + 
  theme_classic() + theme(legend.position="none") +
  ylim(-33.1, -29) + 
  scale_x_continuous(breaks=c(1,3,6,10,20, 40, 60, 100)) 

grid.arrange(p1, p2, nrow=1)

ggsave(plot=p1, width=5, height=3, dpi=300, filename = "dC13_100_paperbirch_L.jpg")
ggsave(plot=p2, width=5, height=3, dpi=300, filename = "dC13_100_paperbirch_R.jpg")

# for white spruce dC13 1-100, use ylim(-33, -28)
# for white spruce dC13 1-20, use ylim(-32.5, -29)
# for paper birch dC13 1-100, use ylim(-33.1, -29)
# for paper birch dC13 1-20, use ylim(-33.1, -29)

#####################################################################
# %N 
#####################################################################
# select data
plant.data <- read.csv("WSU plant isotopes for analysis.csv")

# make distance numeric
plant.data$Distance <- as.numeric(plant.data$Distance)

# by stream 
#plant.data <- plant.data[ which(plant.data$Stream=="Hansen"),]

# by plant 
plant.data <- plant.data[ which(plant.data$Species=="White Spruce"),]
plant.data <- plant.data[ which(plant.data$Species=="Paper Birch"),]

# scale variables
test1 <- scale(plant.data$organic.GWC, scale=FALSE)
plant.data$organic.GWC <- test1[,1]

test2 <- scale(plant.data$mineral.GWC, scale=FALSE)
plant.data$mineral.GWC <- test2[,1]

test2 <- scale(plant.data$mean.mycorrhizal.fungal.C.N, scale=FALSE)
plant.data$mean.mycorrhizal.fungal.C.N <- test2[,1]

test2 <- scale(plant.data$organic.soil.C.N, scale=FALSE)
plant.data$organic.soil.d15N <- test2[,1]

test2 <- scale(plant.data$mineral.soil.d15N, scale=FALSE)
plant.data$mineral.soil.d15N <- test2[,1]

plant.data$decomposing.carcass <- as.factor(plant.data$decomposing.carcass)
plant.data$carcass.deposition <- as.factor(plant.data$carcass.deposition)
plant.data$carcass.load <- as.numeric(plant.data$carcass.load)

# examine effect of distance on dN15 ratio

# for hump, omit transect and mean fungal d15N

# if adding mean mycorrhizal or mean fungal d15N, need to drop rows that do not contain values for these
plant.data <- plant.data %>% drop_na(mean.mycorrhizal.fungal.d15N)
plant.data <- plant.data %>% drop_na(mean.fungal.d15N)

# for white spruce 1-20, have side and distance variables up to fourth power (or just to second) only
# for paper birch, have squared side and distance, plus X.N

# change na. action
options(na.action = "na.fail")
all.parms<-lm(d15N ~ 
                #Side + 
                ln(Distance) + 
                #Side:ln(Distance) +
                #Side:I(ln(Distance)^2) +  
                #Side:I(ln(Distance)^3) + 
                #Side:I(ln(Distance)^4),
                carcass.deposition + 
                decomposing.carcass + 
                carcass.load + 
                X.N + 
                organic.GWC +
                mineral.GWC + 
                mean.mycorrhizal.fungal.d15N +
                #mean.fungal.d15N +
                #mean.fungal.d15N:Side + 
                organic.soil.d15N + 
                mineral.soil.d15N,
              #mean.mycorrhizal.fungal.d15N:Side + 
              #Transect +
              data=plant.data)
results<-dredge(all.parms)
subset(results, delta <2)
subset(results, delta == 0)

# step model
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

# step model with random effedcts
step.model <- lmer(d15N ~ carcass.deposition + ln(Distance) + mineral.soil.d15N +
                     mean.fungal.d15N + organic.soil.d15N + 
                     (1|Stream:Transect), data=plant.data)
summary(step.model)
